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Abstract 

In  this  article  we  present  a  high  resolution  hybrid  central  finite  difference- WENO 
scheme  for  the  solution  of  conservation  laws,  in  particular,  those  related  to  shock- 
turbulence  interaction  problems.  A  sixth  order  central  finite  difference  scheme  is 
conjugated  with  a  fifth  order  weighted  essentially  non-oscillatory  WENO  scheme  in  a 
grid-based  adaptive  way.  High  Order  Multi-Resolution  analysis  is  used  to  detect  the 
high  gradients  regions  of  the  numerical  solution  in  order  to  capture  the  shocks  with 
the  WENO  scheme  while  the  smooth  regions  are  computed  with  the  more  efficient 
and  accurate  central  finite  difference  scheme.  The  application  of  high  order  filtering 
to  mitigate  the  dispersion  error  of  central  finite  difference  schemes  is  also  discussed. 
Numerical  experiment  with  the  ID  compressible  Euler  equations  are  shown. 
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1  Introduction 


Direct  numerical  simulation  of  turbulent  flows  requires  very  accurate  numerical  schemes 
with  the  ability  to  resolve  both  large  and  small  eddies  that  are  often  at  wavelengths  differing 
several  orders  of  magnitude  from  each  other.  High  order  central  finite  difference  (CeFD) 
schemes  are  straightforward,  easy  to  implement,  and  sufficiently  accurate  to  capture  the 
smallest  resolvable  scales  presented  at  these  problems  with  a  small  number  of  points  per 
wavelength.  Many  such  high  order  finite  difference  schemes  had  been  discussed  and  deployed 
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for  solving  linear  advection  problems  such  as  the  Maxwell  equations  in  Electro-Magnetics 
and  wave  equations  in  Acoustics  ([17,4]). 

Another  important  issue  appears  when  dealing  with  the  modeling  of  compressible  flows  via 
the  inviscid  Euler  equations  due  to  the  spontaneous  appearance  of  finite  time  singularities. 
When  applied  to  such  problems,  fixed  stencil  explicit  schemes  like  central  finite  difference 
schemes  produce  oscillations  that  do  not  decrease  in  size  with  grid  refinement.  This  so- 
called  Gibbs  Phenomenon  causes  loss  of  accuracy  and,  if  not  properly  controlled  by  means 
of  artificial  dissipation,  instability  might  follow  as  a  consequence  of  the  contamination  of 
the  solution  by  these  oscillations.  Essentially  Non-Oscillatory  schemes  (ENO)  [14]  have  been 
developed  in  order  to  overcome  the  Gibbs  phenomenon.  ENO  schemes  make  use  of  nonlinear 
weights  based  on  divided  differences  of  the  numerical  solution  in  order  to  bias  the  local 
stencils  when  computing  derivatives,  avoiding  interpolations  across  discontinuities.  Weighted 
Essentially  Non-Oscillatory  schemes  (WENO)  [11,13]  are  an  improvement  over  ENO  schemes 
due  to  their  higher  order  of  accuracy  with  same  stencil  size.  A  convex  combination  of  all 
the  possible  sub-stencils  of  ENO  achieves  optimal  order  of  accuracy  at  the  smooth  parts  of 
the  solution.  Nevertheless,  the  intrinsic  numerical  dissipation  of  WENO  schemes,  although 
necessary  to  properly  capture  shock  waves,  might  seriously  damp  relevant  small  scales. 

A  natural  idea  is  then  to  conjugate  shock  capturing  schemes,  to  be  applied  around  disconti¬ 
nuities,  with  high  order  central  finite  difference  methods  for  the  smooth  parts  of  the  solution. 
In  this  article,  we  study  the  conjugation  of  high  order  central  Finite  Difference  and  WENO 
schemes  when  numerically  solving  systems  of  hyperbolic  conservation  laws.  One  important 
component  is  a  switch  algorithm  to  dynamically  decide, at  any  given  time  step,  which  scheme 
to  turn  on  at  each  grid  point.  The  switch  we  apply  is  based  on  Ami  Harten’s  work  [7,8],  where 
a  high  order  Multi-Resolution  analysis  (MR)  is  performed  at  every  step  of  the  temporal  in¬ 
tegration  process.  The  resulting  adaptive  scheme  ensures  that  fluxes  at  grid  points  around 
discontinuities  will  always  be  computed  by  a  WENO  scheme,  where  as  smooth  tendencies 
will  not  suffer  any  unnecessary  extra  damping  since  they  will  be  treated  by  a  CeFD  scheme. 
Another  advantage  is  that  the  application  of  the  CeFD  scheme  avoids  the  heavy  machinery 
employed  by  the  characteristic-wise  WENO  finite  difference  algorithm  such  as  the  evaluation 
of  the  Jacobian  of  the  fluxes,  characteristics  decomposition  and  recomposition  and  the  global 
Lax-Friedrichs  flux  splitting. 

In  this  work,  we  will  focus  on  the  maximum  order  central  finite  difference  scheme.  Other 
standard  and  dispersion  optimized  central  finite  difference  schemes  can  be  employed  [17] 
but  will  not  be  explored  in  this  study  at  this  time.  Due  to  the  non-dissipative  nature  of 
the  central  finite  difference  scheme,  the  inherent  dispersion  error,  although  small,  generates 
undesirable  oscillations  polluting  the  solution  in  time.  By  applying  high  order  filtering  [16] 
we  are  able  to  mitigate  grid-size  oscillations  that,  otherwise,  would  destroy  the  accuracy. 

The  paper  is  organized  as  follows:  In  Section  2  we  briefly  introduce  the  WENO  finite  differ¬ 
ence  scheme.  Section  3  presents  Harten’s  multi-resolution  algorithm.  A  brief  discussion  on 
the  central  hnite  difference  scheme  is  given  in  Section  4.  The  Hybrid  scheme  is  presented  in 
Section  5.  Numerical  experiments  with  the  proposed  Hybrid  scheme  for  the  one  dimensional 
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Macli  3  shock-entropy  wave  interaction  are  shown  in  Section  6.  The  conclusion  is  given  in 
Section  7. 


2  ENO  and  WENO  schemes 


In  this  section  we  describe  the  ENO  and  WENO  schemes  for  one  dimensional  scalar  conser¬ 
vation  laws  : 

ut  +  f(u)x  =  0.  (1) 

The  conservative  finite  difference  formulation  of  (1)  in  a  uniformly  spaced  grid  is 


duj(t) _ l_  f-  _  f  \ 

dt  Ax  V  'J+\  i-\)  ’ 


(2) 


where  Ax  is  the  grid  size,  Ui(t)  is  the  solution  within  the  stencil  J* 
numerical  flux 


X.  1,1.1 
_  i  2  *+  2  _ 


and  the 


|_1  f{ui—  ri  ui+s)i  (3) 

with  r  and  s,  integer  parameters  defining  the  set  of  values  used  for  the  computation  of  the 
flux  f(u),  satisfies  the  following  conditions: 


•  /  is  a  Lipschitz  continuous  function  in  all  the  arguments; 

•  /  is  consistent  with  the  physical  flux  /,  i.e.,  f(u,u, ... ,u )  =  f{u). 


Thus,  the  Lax-Wendroff  theorem  applies,  and  the  conservative  scheme  (2),  if  it  converges, 
will  converge  to  a  weak  solution. 

The  following  problem  is  related  to  the  computation  of  the  numerical  flux  /  above: 


Given  vt  =  v(xi),  find  v.+  i  =  v(vi-r,  ...,Vi+s)  such  that 


■Ei  (t+i  -  V)  +  0(  Ax’1).  (4) 

In  other  words,  it  is  desired  that  the  flux  difference  above  be  a  fc-th  order  approximation  to 
the  derivative  of  v.  This  can  be  accomplished  by  means  of  an  auxiliary  function  h(x)  which 
can  be  implicitly  defined  as: 


V)  =  ^ 

Since 

^ 

we  take 

v  i  —  h 
*+2  H 

rx-\-  2 

J  *  ww- 

(5) 

Jx 

(h  i  -h.  i  )  , 

V  *+2  1  2  / 

(6) 

i  +0(Axfc), 

(7) 

2 
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and  the  problem  is  solved. 

We  are  now  left  with  the  problem  of  finding  a  fc-th  order  approximation  for  h(x).  We  start  by 
noticing  that,  as  defined  above,  {ty}  are  the  cell  averages  of  h(x)  in  the  intervals  /*.  Therefore, 
the  primitive  function  of  h{x)  can  be  deftned  as  H(x)  =  h(x)dx  and  the  exact  values 
of  H(x)  at  x.  i  is  H[x  i)  =  Ax  ]Tj=-oo  vr  By  interpolating  H(x)  at  these  values,  one  can 

find  a  k  +  1  degree  interpolating  polynomial  P(x).  Thus  the  derivative  p(x)  =  P'(x)  is  the 
desired  /c-th  order  approximation  of  h(x)  we  are  looking  for.  Thus 


X 


i+2 


k- 1 

y  C-rj^i—r+j 


3=0 


h  (xi+i)  +  0{Axk) 


r  —  0, . . . ,  k  —  1. 


(8) 


where  the  coefficients  crj  are  obtained  through  Lagrangian  interpolation  process  (see  [11]). 
They  depend  on  the  order  k  of  approximation  and  also  on  the  left-shift  parameter  r,  but  not 
on  the  values  v*.  Since  r  can  vary  from  0  to  k—  1,  we  have  k  distinct  interpolating  polynomials 
to  choose  from;  all  of  them  yielding  a  fc-th  order  approximation,  once  v  is  smooth  inside  the 
intervals  considered.  The  above  process  is  called  the  reconstruction  step,  for  it  reconstructs 

the  values  of  h  at  the  cell  boundaries  from  the  cell  average  values  h{x)  in  the  interval 

Sr  =  {Uj=o  h-r+i i,r  =  0, . . .  ,k  -  1}  =  { Xi-r ,  ...,xi+s}  with  s  —  k  r  1. 

The  key  idea  of  a  k- th  order  ENO  (Essentially  Non-Oscillatory)  scheme  [14]  is  to  use  the 
smoothest  stencil  Sr  among  the  k  candidates  in  order  to  avoid  Gibbs  oscillations  near  shocks. 
ENO  chooses  the  parameter  r  in  a  way  that  no  discontinuity  is  inside  the  stencil  intervals 
of  the  stencil  Sr.  However,  in  smooth  regions  all  stencils  together  carry  information  for 
an  approximation  of  order  higher  than  r.  WENO  is  an  improvement  over  ENO  for  it  uses 
a  convex  combination  of  all  available  polynomials  for  a  fixed  k,  assigning  essentially  zero 
weights  to  stencils  containing  discontinuities.  This  yields  a  (2 k  —  1)  order  method  at  smooth 
parts  of  the  solution.  The  flux  /.  i  of  the  WENO  method  is  deftned  as 

«+  r) 


k- 1 

r= 0 


(9) 


The  essentially  non-oscillatory  property  is  obtained  by  requiring  that  the  weights  a jr  reflect 
the  relative  smoothness  of  /: 


oor  = 


with 


Cr 

-  0 e  +  ISr)P 


(10) 


where  (3  —  2  and  Cr  are  the  ideal  weights  (see  [2,11])  The  parameter  £  is  set  to  10-10.  A 
good  review  of  the  role  of  the  e  in  the  accuracy  of  the  WENO  scheme  can  be  found  in  a 
recent  paper  [9].  There,  loss  of  accuracy  around  critical  points  in  the  classical  WENO  scheme 
was  also  discussed  and  a  fix  was  proposed,  which  is  also  implemented  in  this  work.  ISj  is  a 
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measure  of  the  smoothness  of  polynomials  on  the  stencils  Sf. 


*■ 


(ii) 


When  the  interpolation  polynomial  prj  in  a  given  stencil  Sj  is  smooth,  the  smoothness  indi¬ 
cator  ISj  is  relatively  much  smaller  than  those  of  stencils  where  the  polynomial  has  discon¬ 
tinuities  in  its  first  r  —  1  derivatives.  Therefore,  discontinuous  stencils  receive  an  essentially 
zero  weight. 

For  the  system  of  Conservation  Laws  such  as  the  Euler  equations  (22),  the  left  and  right 
eigenvectors  and  eigenvalues  of  the  Jacobian  of  the  flux  is  computed  via  the  Roe  Average 
method.  The  global  Lax-Friedrichs  flux  splitting  is  used  to  split  the  flux  into  its  positive 
and  negative  components.  The  resulting  positive  and  negative  fluxes  are  then  projected 
in  the  characteristic  field  via  the  left  eigenvectors.  Each  characteristic  flux  values  is  then 
reconstructed  according  to  the  WENO  algorithm  above  with  one  point  upwinding  bias  at 
each  cell  boundary.  The  numerical  flux  /  i  is  obtained  by  projecting  the  reconstructed 

characteristic  flux  values  back  into  the  conservative  field  via  the  right  eigenvectors  of  the 
Jacobian  of  the  flux.  The  details  of  the  algorithm  can  be  found  in  [11,13]. 


3  Multi-Resolution  Analysis 

It  is  essential  for  the  Hybrid  scheme  proposed  in  this  article  to  measure  the  smoothness  of 
the  solution  in  a  high  order  fashion  such  that  the  appropriate  high  order  algorithm  can  be 
used  at  distinct  grid  points. 

It  is  common  in  the  literature  to  look  at  divided  differences  of  the  solution  in  order  to  indicate 
the  presence  of  discontinuities  at  a  particular  grid  location  when  such  differences  exceed  a 
given  tolerance.  When  performing  the  numerical  calculations  of  Section  6,  the  density  field 
will  be  the  one  used  for  the  analysis,  since  not  only  it  contains  the  discontinuities  due  to  the 
shocks  and  the  rarefaction  waves,  but  also  the  contact  discontinuities,  which  are  present  in 
the  weak  solutions  of  systems  of  conservation  laws  such  as  the  Euler  equations. 

In  this  work,  we  employ  the  Multi-Resolution  (MR) algorithms  by  Harten  [7,8]  to  detect  the 
smooth  and  rough  parts  of  the  solution. 

Consider  a  set  of  nested  dyadic  grids  {Gk ,  0  <  k  <  L},  defined  as: 

Given  an  initial  number  of  the  grid  points  Nq  and  grid  spacing  Axq, 


(12) 
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where  xk  =  jAxk ,  Axk  =  2kAxo,  Nk  =  2  kNo  and  the  cell  averages  fk  of  function  /  at  xk: 


/?  = 


Axi 


y.k 

EJ-1 


f(x)dx, 


(13) 


Let  fkj_ j  be  the  approximation  to  fkj_ l  by  the  unique  polynomial  of  degree  2s  that  interpo¬ 
lates  fk+l,  |/|  <  s  at  xk+l,  where  r  =  2s  +  1  is  the  order  of  approximation  . 

The  approximation  error  (or  multiresolution  coefficients)  clk  =  f^f-i  ~  f'jf-i  i  at  the  k-th  grid 
level  and  grid  point  Xj,  has  the  property  that  if  f(x)  has  p  —  1  continuous  derivatives  and  a 
jump  discontinuity  at  its  p-th  derivative  as  denoted  by  [■],  then 


,k  l  A xpk[f^}  for  p  <  r 

dj~\  ,  ,  •  (14) 

for  p  >  r 

The  approximation  error  dk  measures  how  close  the  data  at  the  finer  grid  ja^-1  j  can  be 
interpolated  by  the  data  at  the  coarser  grid  j.  From  formula  (14)  it  follows  that 

\dkjX  i  ~  2~p\dk\,  where  p  =  min{p,  r},  (15) 

which  implies  that  away  from  discontinuities,  the  MR  coefficients  {dk}  diminish  in  size  with 
the  refinement  of  the  grid  at  smooth  parts  of  the  solution;  close  to  discontinuities,  they 
remain  the  same  size,  independent  of  k. 

The  MR  coefficients  {dk}  were  used  in  [8]  in  two  ways.  First,  finer  grid  data  were  mapped  to 
its  multiresolution  representation  /m  =  {{dk},  •  •  •  ,  {dj},  fL )  to  form  a  multiscale  version  of  a 
particular  scheme,  where  truncation  of  small  quantities  with  respect  to  a  tolerance  parameter 
eMR  decreased  the  number  of  operations  at  fluxes  computation.  Secondly,  the  MR  coefficients 
{dk}  were  also  used  to  generate  a  shock  detection  mechanism  and  an  adaptive  method  was 
designed  where  a  second-order  Lax-Wendroff  scheme  was  locally  switched  to  a  first-order 
accurate  TVD  Roe  scheme,  whenever  {rfl }  was  larger  than  the  tolerance  parameter.  In 
this  article,  we  employ  the  latter  idea  where  a  high  order  central  finite  difference  scheme 
is  switched  to  a  high  order  WENO  scheme  whenever  a  ’’sensor”,  represented  by  the  {dk}, 
detects  discontinuities  at  any  grid  point.  Since  we  will  be  using  the  first  level  k  =  1  only,  we 
shall  drop  the  superscript  1  from  the  dj  from  here  on  unless  noted  otherwise. 

The  adaptive  method  of  Harten  used  only  the  first  level  {dj}  with  a  fixed  MR  order.  However, 
(14)  also  indicates  that  the  variation  of  the  MR  order  can  give  additional  information  on  the 
type  of  the  discontinuity.  For  instance,  consider  the  piecewise  analytical  function 


f(x)  = 


10  +  x3  —l<x<  —0.5 


—0.5  <  x  <  0 


sin(27rx)  0  <  x  <  1 


(16) 
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The  test  function  (16)  has  a  jump  discontinuity  at  x  —  —0.5  and  a  discontinuity  at  its  first 
derivative  at  x  =  0  as  depicted  in  Figure  1  and  is  analyzed  by  the  third  and  the  eighth  orders 
MR.  As  shown  in  Figure  1,  the  quantity  {dj}  at  each  grid  point  decays  exponentially  to  zero 
inside  each  analytical  piece  of  the  function  when  the  order  of  the  MR  is  increased  from 
third  to  eighth.  At  the  discontinuity  x  =  0.5,  {dj}  is  0(1)  and  remained  unchanged  despite 
the  increase  of  the  MR  order.  Similar  behavior  of  the  {dj}  are  exhibited  at  the  derivative 
discontinuity  at  x  =  0  with  a  much  smaller  amplitude.  The  extent  of  the  eighth  order  MR  is 
larger  than  the  third  order  because  more  functional  values  are  needed  to  perform  the  analysis 
at  the  higher  order.  In  general,  NMR  =  2 (nMR  +  1)  number  of  functional  values  are  needed 
for  a  given  MR  order  nMR. 


Fig.  1.  Multi-Resolution  Analysis  of  (Left)  the  test  function  (16)  and  (Right)  the  Lax  problem. 


Also,  in  Figure  1,  the  density  f(x)  =  p(x)  of  the  Lax  shock  tube  problem  with  Riemann  initial 
condition  and  the  MR  coefficients  {dj}  are  shown.  The  location  of  the  shock  at  x  ~  2.5,  the 
contact  discontinuity  at  x  ~  1  and  the  discontinuities  at  the  edges  of  the  rarefaction  wave 
x  ~  —1.8  and  x  ~  —2  are  clearly  delineated  by  the  MR  coefficients  {dj}.  Higher  order 
MR  would  not  make  any  distinguishable  difference  here  as  the  solution  is  a  piecewise  linear 
function. 

Analogy  between  the  wavelet  analysis  and  the  Multi- Resolution  Analysis  is  apparent,  the 
MR  coefficients  { d }}  are  wavelet  coefficients.  Interested  readers  are  referred  to  the  seminal 
book  by  Daubechies  [6]  and  references  contained  therein  for  a  detailed  exposition  of  Multi- 
Resolution  Analysis  in  the  context  of  wavelet  analysis. 


4  Central  Finite  Differences  Schemes 


In  the  following  discussion,  the  computational  grid  is  restricted  to  a  uniformly  spaced  grid, 
i.e.  {xj  =  j  Ax,  j  =  0, . . . ,  N},  where  Ax  is  the  grid  spacing.  For  a  given  data  {fj  =  f(xj),j  = 
0, . . .  ,N},  the  first  derivative  of  the  function  j-  f(xj )  can  be  approximated  to  order  n  using 
2n  +  1  grid  points  { Xj-n , . .  .  ,£?■+„}  and  its  corresponding  functional  values  {/)_„, . .  . ,  fj+n} 
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as 


Txf(Xi)  =  7TX^-  a ih+i' 


(17) 


3=~n 


where  w3  are  the  Lagrangian  weights  of  the  first  derivative  and  here  we  assumed  N  >  2n  +  l. 


One  of  the  good  properties  of  a  central  finite  difference  scheme  is  that  it  is  non-dissipative  in 
nature;  however,  dispersion  is  also  present,  requiring  the  use  of  high  order  smoothing.  The 
non-dissipative  property  is  guaranteed  by  requiring  w-j  =  —Wj.  To  examine  its  dispersive 
property,  the  Fourier  transform  is  applied  to  (17),  yielding  the  effective  wavenumber  k  of  the 
scheme,  which  can  be  expressed  in  terms  of  the  exact  wavenumber  k  in  the  form: 

n 

kAx  —  2  y~]  Wj  sin(jkAx).  (18) 

3= 1 

The  difference  between  the  effective  wavenumber  k  and  the  exact  wavenumber  k  is  the 
dispersion  error  E(kAx)  =  \kAx  —  kAx\.  The  higher  the  order  n  of  the  scheme,  the  better 
the  scheme  can  resolve  higher  wavenumbers.  However,  the  error  analysis  also  indicates  that 
the  wavenumber  k  =  can  never  be  resolved  regardless  of  the  order  of  the  scheme  (see 
Figure  2).  In  other  words,  without  artificial  dissipation,  one  can  expect  high  wavenumbers 
(frequencies)  error  in  the  solution  of  a  nonlinear  hyperbolic  PDE. 


To  mitigate  dispersion  error  on  the  solution,  we  propose  to  apply  high  order  filtering.  By 
this,  we  mean  that  the  large  scales  of  the  function  are  preserved  up  to  the  order  of  the  central 
finite  difference  scheme,  while  the  small  scales  usually  associated  with  noise  and  error  will 
be  attenuated  smoothly  to  zero. 

For  a  given  function  /,  discretized  on  a  uniformly  spaced  grid,  the  filtered  function  /  at  the 
grid  point  x*  can  be  expressed  as 

n 

ft  =  J2  ajfc+j,  (19) 

j=-n 

where  aq  are  the  filtering  weights  of  order  n.  The  filtering  weights  satisfy  the  symmetry 
property  (\-3  =  oq,  ensuring  no  dispersion.  The  a  are  chosen  in  such  a  way  that  the  Erst 
n  moments  of  the  hltered  function  match  exactly  the  first  n  monomials  {1,  x,x2, . .  .  ,xn } 
ensuring  that  the  approximation  order  of  the  hltered  function  is  kept  high.  In  addition  to 
that,  the  a  are  also  required  to  satisfy  the  condition  Yl]=-n  =  0  so  that  oscillations 

at  high  wavenumbers  k  near  are  attenuated  to  zero.  By  applying  the  Fourier  transform 
to  (19),  the  physical  filter  can  be  expressed  as  a  filter  cr(kAx)  in  the  frequency  space, 

n 

a(kAx)  =  «o  +  2  cos(jkAx),  (20) 

3= 1 

with  cr(0)  =  1  and  cr(7r)  =  0  for  normalization.  The  higher  the  order,  the  hatter  is  the 
hlter  near  the  wavenumber  k  =  0  (see  Figure  2).  Therefore,  more  low  wavenumbers  remain 
unaltered  and  less  higher  wavenumbers  get  attenuated  to  zero.  Some  of  these  filtering  weights 
a  can  be  found  in  [16]. 
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Fig.  2.  The  dispersion  relation  (left)  and  filter  a(kAx)  (right)  of  Central  finite  difference  scheme  of 
orders  n  =  4,  6,  8, 10. 

5  High  Order  Hybrid  Central  -  WENO  Finite  Difference  scheme 


In  this  section,  we  introduce  the  High  Order  Hybrid  Central  Finite  Difference- WENO  finite 
difference  scheme  (Hybrid)  for  the  Euler  equations  and  discuss  related  issues  such  as  the 
dispersion  error  associated  with  it.  In  this  particular  context,  we  define  the  Hybrid  scheme 
as  a  grid-based  adaptive  method  in  which  :w  the  choice  of  the  numerical  scheme  is  determined 
by  the  smoothness  of  the  solution  at  each  grid  point. 

Like  many  adaptive  methods  in  the  literature,  the  smoothness  of  the  solution  is  measured  by 
some  criteria,  such  as  the  gradient  of  the  solution  [1,12],  or  by  the  Multi-Resolution  procedure 
[3,5].  A  finite  difference  scheme  such  as  a  Compact  scheme  [1]  or  an  optimized  central  finite 
difference  scheme  [10]  can  be  used  at  those  grid  points  where  the  solution  is  flagged  as  smooth 
in  lieu  of  the  standard  high  order  WENO  scheme.  By  avoiding  the  computationally  expensive 
characteristics  decomposition/recomposition  and  Jacobian  evaluations  of  the  WENO  finite 
difference  scheme,  it  is  anticipated  that  the  Hybrid  scheme  will  provide  greater  efficiency 
and  resolution,  in  particular,  for  large  scale  simulations  of  turbulence  [10]. 

Here  are  some  details  of  the  implementation  of  the  Hybrid  scheme: 

•  The  Hybrid  scheme  consists  of  a  sixth  order  Central  Finite  Difference  scheme  (CeFD6)  and 
a  fifth  order  WENO  finite  difference  scheme  (WEN05)  for  the  smooth  and  discontinuous 
parts  of  the  solution  respectively. 

•  The  third  order  three  stage  Runge-Kutta  TVD  scheme  (RK3)  is  employed  for  the  temporal 
evolution  of  the  resulting  ODE  with  CFL=0.4  (see  [11]). 

•  The  Multi-Resolution  Analysis  is  applied  only  once  at  the  beginning  of  the  Runge-Kutta 
time  stepping  scheme.  A  grid  point  is  flagged  as  non-smooth  when  the  MR  coefficient 
\di\  >  eMR,  a  user  defined  MR  tolerance  and  the  WENO  flux  will  be  employed  there. 


Flag* 


1, 

0, 


\di\  >  e 


MR 


otherwise 


(21) 
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•  Once  the  flags  are  set,  a  number  of  neighboring  points  around  each  flagged  point  ay, 
depending  on  the  number  of  the  ghostpoints  needed  for  a  given  order  of  the  CeFD  scheme 
and  the  WENO  scheme,  are  also  flagged  to  1. 

In  particular,  if  Nc  and  Nw  are  the  orders  of  the  CeFD  scheme  and  the  WENO  scheme 
respectively,  the  number  of  ghostpoints  required  by  the  CeFD  scheme  and  the  WENO 
scheme  are  |iVc  and  ^(Nw  +  1)  respectively.  At  any  given  point,  say  ay,  flagged  as  non¬ 
smooth,  its  r  =  max  (^Nc,  |  (Ay,  +  1))  neighboring  points  {ay_r, . . . ,  ay, . . . ,  ay+r}  will  also 
be  designated  as  non-smooth,  that  is,  {Flag^  =  1  ,j  =  i  —  r, ...  ,i  +  r}.  This  procedure 
avoids  computing  the  derivative  of  the  solution  by  the  CeFD  scheme  using  non-smooth 
functional  values.  Furthermore,  the  same  WENO  flag  will  be  used  at  all  the  Runge-Kutta 
stages  and  will  be  updated  at  the  next  time  step. 

•  Since  the  CeFD  scheme  is  numerically  more  efficient  when  compared  to  the  WENO  scheme, 
we  opt  for  initially  computing  the  flux  over  the  complete  flow  field  with  the  CeFD  scheme. 
The  WENO  scheme  is  then  used  to  overwrite  the  flux  at  those  grid  points  designated  as 
non-smooth  by  the  WENO  flag. 

•  The  class  of  problems  studied  are  restricted  to  those  where  the  boundary  conditions  do 
not  present  any  complications  with  the  ghostpoints,  for  instance,  periodic  or  freestream 
boundary  conditions.  We  shall  use  as  many  ghostpoints  as  required  for  the  given  order  of 
the  CeFD  scheme,  the  WENO  scheme  and  the  Multi- Resolution  Analysis. 


6  Numerical  Experiments 


Let  us  now  consider  the  one  dimensional  Euler  Equations  for  gas  dynamics  in  strong  conser¬ 
vation  form: 


dW  dF 
dt  dx 


(22) 


W  =  \p,  m,  E]T  is  the  vector  of  conservative  variables  and  the  flux  F  is  given  by 


F 


171 

m, - b  P, 

P 


( E  +  P ) 


(23) 


Here  p,m,  and  E  are,  respectively,  the  density,  mass  flux  and  total  energy  per  unit  volume. 
It  is  coupled  with  the  equation  of  state  for  ideal  gas,  P  —  (7  —  1)(E  —  \m2 / p)  with  7  =  1.4. 


Consider  the  one  dimensional  Mach  3  shock-entropy  wave  interaction,  specified  by  the  fol¬ 
lowing  initial  conditions: 


f  (3.857143,2.629369,10.33333)  x  <  -4 
(P,  «,  P)  =  ,  (24) 

I  (1  +  esin(kx),  0, 1)  x  >  — 4 


where  x  £  [—5,  5]  ,  e  =  0.2  and  k  —  5.  The  solution  of  this  problem  consists  of  not  only  large 
scales  waves  but  also  shocklets  and  fine  scales  structures  which  are  located  immediately 
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behind  the  shock.  We  computed  the  solution  at  time  t  =  1.8  with  N  =  400  uniformly  spaced 
grid  points  using  both  the  standard  fifth  order  WENO  finite  difference  scheme  (WEN05) 
and  the  Hybrid  sixth  order  central  finite  difference-fifth  order  WENO  finite  difference  scheme 
(Hybrid).  The  parameters  used  in  the  Multi-Resolution  Analysis  for  the  Hybrid  scheme  are 
order  nMR  =  8,  tolerance  eMR  =  5  x  10~3  using  the  density  p  as  the  test  function.  In  the  left 
figure  of  Figure  3,  the  solid  black  line  and  red  square  symbols  depict  the  density  computed 
with  the  WEN05  scheme  and  the  Hybrid  scheme  respectively.  The  red  dash  line  encloses 
the  grid  points  where  the  solution  is  deemed  insufficiently  smooth  by  the  MR  Analysis.  The 
WEN05  scheme  is  employed  at  these  grid  points  in  order  to  capture  the  main  shock  and 
shocklets  and  avoid  oscillations.  The  CeFD6  scheme  is  used  for  the  rest  of  the  grid  points.  At 
this  scale,  the  results  of  the  WEN05  scheme  and  the  Hybrid  scheme  are  in  good  agreement 
with  each  other. 


Fig.  3.  (Left)  :  The  density  solution  of  the  Mach  3  shock-entropy  wave  interaction  with  N  =  400  at 
time  t  =  1.8.  (Middle)  :  Zooming  of  the  inflow  region  and  dispersion  oscillations.  (Right)  :  Shows 
the  effectiveness  of  the  filtering  in  controlling  the  dispersion  error.  The  solid  line  and  the  square 
symbols  are  solutions  as  computed  by  the  WEN05  scheme  and  the  Hybrid  scheme  respectively. 

As  mentioned  earlier,  wave  dispersion  is  an  issue  when  using  a  non-dissipative  central  finite 
difference  scheme.  When  zooming  at  the  inffow  region  x  G  [—5,  —3]  (middle  figure  of  Figure 
3),  one  observes  the  oscillations  caused  by  dispersion.  As  proposed  in  last  section,  we  apply 
high  order  filtering  in  order  to  mitigate  the  dispersion  error.  The  order  of  the  filter  nF  used 
in  this  study  is  always  2  orders  higher  than  the  order  of  the  central  finite  difference  scheme 
nCFD,  i.e.  nF  =  nCFD  + 2.  Furthermore,  the  filtering  operation  will  be  performed  at  the  end  of 
the  Runge-Kutta  time  step.  Filtering  in  the  intermediate  stages  of  the  Runge-Kutta  scheme 
is  also  possible  if  deemed  necessary. 

The  density  solution  of  the  shock-entropy  wave  interaction  as  computed  by  the  filtered  version 
of  the  Hybrid  scheme  and  the  zooming  near  the  inflow  region  are  shown  at  the  right  figure 
of  Figure  3.  It  can  be  easily  observed  that  the  eighth  order  filtering  improves  the  quality  of 
the  solution.  The  filter  effectively  controls  the  dispersion  error  and  localizes  the  oscillations 
close  to  the  steep  gradients. 

We  now  illustrate  the  algorithm  with  a  longer  time  simulation  by  computing  the  solution  up 
to  time  t  —  5  in  a  larger  physical  domain  x  G  [—5,  —15]  and  N  =  800  with  both  the  WEN05 
scheme  and  the  Hybrid  scheme  with  MR  tolerance  eMR  =  10-2.  Density  plot,  zooming  of 
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the  high  frequencies  solution  behind  the  main  shock  and  the  MR  coefficients  {dj}  are  shown 
at  time  t  =  5  in  Figure  4.  The  main  shock  and  the  shocklets  are  well  tracked  and  one  sees 
that  the  small  scales  behind  the  main  shock  are  better  resolved  with  the  Hybrid  scheme 
than  with  the  WEN05  scheme.  Notice  that  given  an  appropriate  threshold  for  the  switch  of 
the  schemes,  the  coverage  of  the  WEN05  scheme  remains  relatively  quite  small,  indicating 
that  the  solution  is  mainly  computed  by  the  more  efficient  CeFD6  scheme.  Furthermore,  the 
decay  of  the  MR  coefficients  {dj}  near  the  inflow  region  also  indicates  the  fast  decay  of  the 
dispersion  error. 


Fig.  4.  (Left)  :  The  density  solution  of  the  Mach  3  shock-entropy  wave  interaction  with  N  =  800 
at  t  =  5.  (Middle)  :  Zooming  of  the  high  frequency  waves  behind  the  main  shock.  (Right)  :  The 
multiresolution  coefficients  dj  at  t  =  5. 

The  Hybrid  scheme  has  the  advantage  of  being  considerably  faster  than  the  WEN05  scheme. 
Even  though  timing  of  the  one  dimensional  system  of  equations  is  not  representative  as  timing 
of  a  large  scale  simulation  of  a  multi-dimensional  system  of  PDEs,  we  provide  a  table  listing 
the  average  CPU  usage  in  seconds  per  each  Runge-Kutta  step  of  the  WEN05  scheme  and 
the  Hybrid  scheme  just  to  illustrate  the  potential  gain  in  computational  speed.  From  table  I, 
the  gain  in  speed  can  be  as  much  as  a  factor  of  2  for  a  number  of  grid  points  N  =  800, 1000 
and  up  to  as  much  as  a  factor  of  5  for  larger  N  =  2000.  It  is  obvious  that  the  speedup 
is  a  function  of  many  factors  such  as  the  physical  problem,  the  structure  of  the  solution, 
the  choices  of  the  parameters,  the  computational  platform  and  programming  efforts.  A  more 
detailed  study  on  computational  time  of  the  Hybrid  scheme  in  comparison  with  the  standard 
WENO  scheme  will  be  given  in  a  future  paper. 


N 

WEN05 

Hybrid 

Speedup 

800 

3.8e-3 

1.6e-3 

2.4 

1000 

5.1e-3 

2.1e-3 

2.4 

1500 

1.3&-2 

3.2e-3 

4.1 

2000 

2.6e-2 

4.7e-3 

5.5 

Table  I 

Average  CPU  time  (seconds)  per  Runge-Kutta  step  as  a  function  of  N.  The  last  column  shows  the 
speedup  factor  by  the  Hybrid  scheme  over  the  WEN05  scheme. 
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7  Conclusion 


For  the  high  resolution  solution  of  the  nonlinear  conservation  laws,  we  have  conjugated  an 
efficient  and  accurate  high  order  central  finite  difference  scheme  and  the  high  order  shock- 
capturing  WENO  finite  difference  scheme  together  to  form  a  hybrid  scheme.  The  spatial 
and  temporal  hybridization  of  these  two  schemes  are  accomplished  via  the  high  order  Multi- 
Resolution  analysis.  The  numerical  error  due  to  the  dispersive  nature  of  the  central  finite 
difference  scheme  is  controlled  and  localized  by  the  employment  of  the  high  order  finite 
difference  filter,  which  maintains  the  formal  high  order  of  accuracy  of  the  hybrid  scheme. 
The  main  advantages  of  the  hybrid  scheme  are  the  potential  large  reduction  of  CPU  time 
needed  to  perform  large  scale  numerical  experiments  and  improvement  in  overall  accuracy 
over  the  classical  high  order  WENO  finite  difference  method.  The  reduction  of  the  CPU  time 
usage  is  mainly  due  to  removal  of  calculation  of  the  Jacobian  of  the  flux  of  the  PDE,  the 
characteristic  decomposition  needed  to  project  the  flux  from  the  conservative  held  into  the 
characteristic  held,  where  the  WENO  reconstruction  of  the  numerical  hux  takes  place,  and 
back,  in  the  smooth  regions  of  the  how  held.  The  improvement  of  the  overall  resolution  is 
due  to  the  non-dissipative  nature  of  the  central  finite  difference  scheme  that  yields  a  better 
representation  of  the  solution  over  the  upwinding  tendency  of  the  WENO  reconstruction 
in  order  to  avoid  using  stencils  with  potential  high  gradients.  The  success  of  the  hybrid 
scheme  depends  very  much  on  the  properly  tuning  of  the  parameters  associated  with  the 
Multi-Resolution  analysis  and  is  a  subject  of  the  research  in  progress. 
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